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Abstract 

Coherent energy transfer in pigment-protein complexes has been studied by mapping the quantum network 
to a kinetic network. This gives an analytic way to find parameter values for optimal transfer efficiency. In the 
case of the Fenna-Matthews-Olson (FMO) complex, the comparison of quantum and kinetic network evolution 
shows that dephasing-assisted energy transfer is driven by the two-site coherent interaction, and not system-wide 
coherence. Using the Schur complement, we find a new kinetic network that gives a closer approximation to the 
quantum network by including all multi-site coherence contributions. Our new network approximation can be 
expanded as a series with contributions representing different numbers of coherently interacting sites. 

For both kinetic networks we study the system relaxation time, the time it takes for the excitation to spread 
throughout the complex. We make mathematically rigorous estimates of the relaxation time when comparing 
kinetic and quantum network. Numerical simulations comparing the coherent model and the two kinetic network 
models, confirm our bounds, and show that the relative error of the new kinetic network approximation is several 
orders of magnitude smaller. 

Keywords: exciton transfer, quantum efficiency, kinetic networks, FMO, coherent energy transfer, quantum 
networks, Schur complement. 
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1 Introduction 

Since coherent energy transfer in the Fenna-Matthews-Olson complex (FMO) has been observed [6l [9j [T3] , extensive 
experimental and theoretical research has been dedicated to studying coherent resonant transfer [5] and the coherent 
pigment-protein interaction[12, 8J. In particular, numerical solutions of simple models have shown that dephasing - 
the destruction of the coherences - at an intermediate rate helps to increase the energy transfer efficiency |1Q[ [TT] . 
This has been called dephasing- or environment-assisted energy transfer, and is analogous to a critically damped 
oscillator. The dephasing corresponds to damping and causes the exciton to relax to an equal distribution for every 
pigment site instead of staying localized due to the energy mismatch between the sites. 

The models are based on two assumptions. First, only a single exciton is present, it is located at any of the seven 
pigments. The pigment exciton energy, and the pigment dipole-dipole interaction [H [T] then lead to an oscillatory 
evolution of the system. And second, the site-environment interactions are assumed to be purely Markovian without 
any temporal or spatial correlations. The environment interactions are dephasing, recombination and trapping. 
Dephasing destroys the site coherences without destroying the exciton itself, and phonon recombination or photon 
re-emission lead to loss of the exciton to the environment. Trapping is the transfer of the exciton to the reaction 
center, where the electronic energy is converted to chemical energy, in FMO it occurs at pigment 3. The transfer 
efficiency is the probability that an exciton starting at site 1 or site 6 reaches the reaction center. For a general 
system with n pigments, we convert the master equation of the coherent model into vector form 

P = M P 

where p <E W 1 is the density matrix in vector form and M is a real n 2 x n 2 -matrix. Two procedures to find M are 
presented in 12.21 and 13.51 
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Figure 1: The efficiency of energy transfer in the FMO monomer. Model parameters are described in 17. 31 

To study population transfer channels and conditions for optimal transfer, a mapping to kinetic networks has 
been proposed [3j[7]. A kinetic network is a system where the exciton jumps incoherently between sites according 
to some fixed rates, i.e. a continuous-time Markov process. In its simplest version this approximation only takes 
into account the coherent interaction between pairs of sites to derive the transfer rate between them. If the two 
sites interact with strength V, have an energy separation E, and both sites experience dephasing at rate 7 and 
population loss at rate k then the rate is 

This rate is maximized for the intermediate dephasing rate 7 = E — k so the phenomenon of dephasing-assisted 
transfer is maintained in this approximation. For a system with n sites, these rates constitute the off-diagonals of 
a n x n rate matrix TVo, and the system populations evolve according to 

P = N p 

where p G R n is the time-dependent population vector. Figure [1] displays the transfer efficiency with models M and 
No for different 7, the dephasing-assisted regime clearly shows as a peak around 7 « 170cm -1 . At the peak the 
population evolution of M is well approximated by that of No, therefore dephasing-assisted energy transfer can be 
explained by the relatively simple coherent dynamic between pairs of sites that enters the rate /i and the influence 
of system-wide coherence is small. 

To extract the limit of good approximation we introduce scaling variables, T which is proportional to the energy 
separations, dephasing and population loss rates, and O which is proportional to the site interactions. We will show 
that the approximation of No to M becomes good as Or -1 approaches 0. We generalize the procedure of finding a 
kinetic network approximation in a mathematically appealing way using block matrices. We find a kinetic network 
matrix TV that follows the evolution of M much closer - it is over three orders of magnitude more precise than the 
network No as shown in Figure [TJ Further, it can be expanded in Or -1 as 
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where No is the approximation described above, and the Nk are rate corrections due to coherent interactions via 
k intermediate sites. The expansion terms become smaller for increasing k, N k oc 9 • (6r~ 1 ) fc . By stopping the 
expansion at a finite k kinetic networks approximation of varying accuracy can be formed allowing the study of 
coherent interaction at different "scales" or number of involved sites. We restrict our further investigation to the 
dominant contribution No and the entire sum N. 

In our exact bounds we study the system with all population-loss mechanisms removed. Due to dephasing the 
exciton spreads throughout the system at the exciton relaxation time r and all populations become equal. The 
difference At between relaxation times of M and N or No gives a simple measure of how good the kinetic networks 
approximate the quantum network. As 0r _1 becomes small the kinetic networks approach the quantum network 
and At becomes small as well. 

We define t and At as follows, using the Euclidean norm \\p\\ 2 — y/Y^ii=i Pi t° compare population vectors. 

Definition 1. 

1 . The map T : W 1 —> R" is the restriction of density vectors p to population vectors p, and consequently 
gives the embedding of population vector space in density vector space. In particular, if the first n components 
of p represent the site populations, then T — (1„, nx („2_ n )). 

2. The maximum relaxation time is 

F°° 1 

r = max / e Nt po alt 

pa Jo n 



and the corresponding minimal relaxation rate is 



1 

M = - 

T 



The maximum deviation of relaxation time between the quantum network M and the kinetic network N is 



At = max 

Po 



/>oo 

/ (Te Mt T^-e Nt )p dt 
Jo 



4. Define t , fio an d Ato in the same way, replacing N with N . 

For our bounds we require that every site experiences dephasing. Further, the network has to be connected, 
meaning that any two sites can exchange populations -directly or indirectly- such that the relaxed state will have 
equal population everywhere. And finally we also require our site interactions to be real - but it is clear from our 
proofs that the generalization to complex interactions could be treated in a similar manner. 

Our first results shows how fast the relaxation time of the two kinetic networks No and N approximate that of 
the quantum network M as Or -1 gets small. 

Theorem 2. There are scaling invariant constants ki and ki, such that for QT^ 1 small enough we have the 
following bounds: 

1. The relative difference of relaxation time between quantum evolution M and kinetic evolution No is bounded 
by 

AT ,rci = Ato/to < fcier- 1 . 

2. The relative difference of relaxation time between quantum evolution M and kinetic evolution N is bounded 
by 

Ar re ] = At/t < fc 2 e 2 r~ 2 . 

This Theorem follows from Theorem [5] and Corollary [7] in Section [5l 
We also find the following exponential bounds on the time dependence. 

Theorem 3. There are scaling invariant constants k%, ki and k§, such that for any initial population distribution 
Po we have the following bounds, as long as QT^ 1 is small enough: 

1. For all times t > 

\\Te Mt Tip - e Not po\\ 2 < k 3 e-^ 2 ■ QT' 1 . 

2. For all times t > 

\\Te Mt T^p Q - e Nt ph\\ 2 < k^' 2 ■ 9 2 r- 2 (l + fcslogGr- 1 ) . 

This Theorem follows from Theorem|5]and Corollary[Tn]in Section[5] We expect that more sophisticated methods 
might yield the same bound without the 6 2 r _2 log6r _1 term. 
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2 The quantum network 



We first introduce the Master equation for the coherent model. Then we reformulate the equation in vector form and 
combine the entire dynamic in the real n 2 x n 2 -matrix M. We describe the general structure of M as a preparation 
to the next section, where we generate kinetic networks from parts of M. 

2.1 Master equation 

We consider the same quantum mechanical system studied in with n sites carrying a single excitation which is 
equivalent to a system with n states/levels. The site energies are £ K so the energy operator is 

n 

H = Y / E k \k)(k\. 

k=l 

The site k couples to site I with interaction strength Vki 6 C so the interaction operator is 

V = J2Vu\k)(l\ 

kjtl 

where Vki = Vik- Site trapping, re-emission and recombination can be incorporated by an anti-hermitian operator 
A. Let Kk be the combined rate of exciton loss at site k due to these effects, then A is defined as 

n 

A=^^K k \k)(k\. 
k=l 

Finally, every site is also under the influence of dephasing at rate 7^ > incorporated in the Lindbladian superop- 
erator 

n 

C(p)=Y,LkpLl--{p,LlL k } 

k=l 

with Lk = y/jk\i)(i\. Setting h = 1, the single exciton manifold of the quantum network is described by the master 
equation 

p=-i[H + V, P ]-i{A,p} + £(p) (2) 

where square and curly brackets represent commutator and anti-commutator respectively. 

For now we set A = 0, ignoring exciton depleting processes as explained above. We will mention how to include 
them in the kinetic network approximations later on. Our approximation becomes exact in the limit where the 
energy difference between sites is large, the dephasing is large and the interactions are small. To be specific, we 
introduce scaling parameters T and O and consider the limit 9r _1 — > 0. Energies and dephasing scale like T and 
interactions scale like 

Ek^r, 
ik oc r , 
Vki oc e . 

With these assumptions the master equation turns into 

p=-i[TH + eV,p}+TC(p). (3) 
Because this equation is linear in p it can be converted into vector form 

'p = Mp 

where p G M™ is the density matrix in vector form and M is a real n 2 x n 2 -matrix. Two procedures to find M are 
presented in 12.21 and 13.51 
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2.2 Converting to vector equation 

We rewrite the master equation ([3]). skipping the scaling factors and T, it is easy to reintroduce them at a later 
point 

p = -i[H + V,p]+C(j>). (4) 
Our first goal is to convert this into the differential equation 

'p = Mp 

for density "vectors" p € K" 2 . Notice that because p — p^ the space of density matrix has n 2 real dimensions, so we 
are not using any information when mapping p to p. 
We use the following conversion: 

1 . The first n entries of the density vector are the populations - the real diagonal entries of p. 

2. For the entries n + 1 to n 2 we alternate between real and imaginary parts of the coherences - the off-diagonal 
entries of p - starting with entry p^i where k — 1 and I — 2 continuing by increasing I until I = n, then moving 
to the entry /?23- We multiply all these entries with \^2, a normalization factor useful to achieve simpler 
expressions later on. 

In terms of index equations this is: 

1 . For k = 1 . . . n 



2. For k, I € {1, . . . , n} with k < I 



Pk = Pkk ■ 

Pn+2n(k-l)-k(k+l)+2l-l = ^2 Re Pkl 
Pn+2n(fc-l)-fc(fc+l)+22 = V2Im Pkl ■ 

Other mappings will yield the same kinetic networks, as long as they allow for an easy separation of population 
and coherence space. 

While somewhat tedious, it is now relatively straightforward to find the matrix M such that 

'p = Mp. 

To find the rows k = 1 ... n we write out the diagonal components of the RHS of (Q}, and to find rows k = n+1, . . . , n 2 
we write out the off-diagonals of the RHS of (Q}. We follow this procedure explicitly for the case n = 3 in Appendix Rl 
From there it is obvious how the procedure generalizes to larger n. Here we will only present the final form. 

2.3 The coherent evolution matrix M 

2 

For simple notation and to simply extract the kinetic networks we split up the density vector space W 1 . Let P = M 
be the space of populations and let C = M" ~™ be the space of coherences. We can then write density vectors as 



with p G P and c£C. With this splitting the matrix M describing the quantum network looks like 

'0 -a tN 



M 



where a : P — > C and b : C — > C are real matrices (so = a T , but we'll keep the more general notation for later). 
Notice that the populations do not affect each other directly, but only via the coherences. 

Matrix a describes how populations couple to coherences, its entries are real and imaginary parts of Vki, naturally, 
site k will only couple to coherences kl with I ^ k, thus of the (n 2 — n) entries in the fc-th column of a only 2(n— 1) 
are nonzero. Matrix b describes how coherences couple to other coherences, if considered as a block matrix with 
2 x 2-blocks the diagonal block for the coherence between site k and site I is of the form 



~lkl —Efel 

Eh — 7fe( 



(5) 
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where 



1m = ^(ik+li) 
Em = Ek — Ei . 

The off-diagonal blocks consist of real and imaginary parts of Vu ■ 

From the form of M, when ignoring the off-diagonal blocks of b, we see that the site k couples to the site I via 
the coupling strength Vu, then some mixture of ju and Eki and then again via the coupling strength Vu- This 
reminds us of the rates of the form fi = -x^gs described in (fl) that make up the matrix N . We will make this 
intuition precise is the next subsections. 



3 Kinetic networks 

In this section we show how the kinetic network N emerges naturally out of the study of the resolvent (z — M) _1 . 
We expand N in powers of Qr _1 , giving the series 

oo 

k=a 

with the leading order contribution being N . For some steps involving matrix calculations we only give a simplified 
version. However, in Appendix (|A"|) we follow the procedure described below, giving the full expressions in the case 
n = 3. 



3.1 Extracting the kinetic network N 

To extract kinetic networks from M we consider its resolvent (z — M) . Remember that for any holomorphic 
function / we have 

f(M) = ±-j)f(z)(z-M)- 1 dz. 

Therefore, if one can bound the resolvent appropriately, one can also bound the evolution operator e Mt and other 
related quantities. Because we only care about approximating the population dynamics we restrict our view to the 
population block of the resolvent of M . The Banchiewicz formula |2] gives the inverse of a 2 x 2-block matrix. 
The first block of the inverse - in our case the population block - is called the Schur complement, and due to its 
basic nature has many applications in applied mathematics, statistics and physics |14) . Here we use it to "pull" the 
coherence dynamic back into population space. Only writing the Schur complement and skipping the other blocks 
of the resolvent we have 

(z-M)- 1 = (^ z - a Hb- z^a)- 1 A 

Remember the operator T, the restriction to population space. With our choice of density vector basis it has the 
form 

T = (l n nx („2_ n )J . 

The difference of evolution for initial conditions p = Pq°^ = ^Po (zero coherences) between quantum network M 
and kinetic network N is thus 

( Te Mt r t _ e m\ p Q = _L I e zt n z _ a ]r h _ z) -i a) -i _ {z _ N) n dz . 
v 2ni J 

For a good approximation we require 

(z-o)(b- z^a)' 1 w (z- N)K (6) 

At this point it is a small step to drop the second z on the LHS, in which case the formula becomes equality if we 

set 

N = a%~ 1 a. 
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To see intuitively that this approximation is good, consider the following. Matrix b contains terms proportional to 
r on its diagonal and terms proportional to on its off-diagonal, matrix a is proportional to 0, therefore 

n oc e 2 ^ 1 

when 8r -1 becomes small. For values of z that are smaller than eigenvalues of b the approximation © is good 
because (b — z)^ 1 m for values of z larger than eigenvalues of b is good, because then z is much larger than the 
eigenvalues of N, and so both sides of ((6j) are approximately z _1 . This basic insight is what drives our bounds in 
Section M 



3.2 Expanding N 



As mentioned in 12. 31 b consists of 2 x 2-blocks proportional to T on the diagonal and 2 x 2-blocks proportional to 
on the off-diagonal. We separate this contributions defining 

b = b + v 

where bo oc T and v oc is the block-diagonal and block-off-diagonal of b respectively. If 0r _1 is small enough and 
if 6o is invertible we can expand 

oo 

. k 



This leads to the expansion 



with 



b- 1 =^ (-vb^) 

k=0 

oo 

N = s b' 1 a = J2 Nk 



fc=0 



N k = a% 1 {-ub Q - 1 ) k a. (7) 



When using explicit forms of a, bo and v one can see that the rates in N k consist of corrections due to interactions 
via k intermediates. Roughly speaking, every of the [k + 1) sites along the chain contributes a factor of 0, every of 
the k coherences (links) contributes a factor of thus N k scales like fc+1 r~ fe . 



3.3 The network N 

We now present the explicit form of 



iV = a f 6 1 a 



the dominant contribution to TV. We only show the crucial parts of the calculations that should make clear how to 
get the result for general n. 

Notice that, from 13.21 and it follows that b is a (n 2 — n) x (n 2 — n) matrix with the only nonzero entries 
being 2x2 blocks 

Eki — 7fci 

along the diagonal. With the unitary transformation 

we can diagonalize these 2x2 blocks. Hence, the entire matrix bo can be diagonalized by applying the transformation 

U = l(„2-„)/ 2 ® Uo (8) 

and 

b = U^boU = diag(ai2, ai2, ai3, a i3 , a n -i,n) (9) 

with 

a-ki = -7fci + iE ki 
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C/ f o 



(10) 



where diag denotes a diagonal matrix with given diagonal entries. In fact, U also helps to simplify a, consider the 

case n = 3 

/Via -Via \ 

Via -Via _ 

V13 -V13 

V13 _ -Via 

V23 —V23 

V v 23 -V23J 

and the same happens for v = U^uU (derivation in Appendix [A"]) . Notice that both 60 and a are complex matrices, 
still, we can use the transformed matrices a, 60, and v when finding explicit expressions for the real matrices 
because U cancels out. For example 



In the case n = 3 we get 



with 



N = alb^a 

= alUtfb^Utfa. 



-Mia — Mis M12 Mis 

A*12 — Ml2 _ ^23 M23 

Ml3 A*23 — Ml3 - A*23; 



(11) 



Hki 



2|Vfc;| 2 7fci 



The following simplified calculation illustrates how the rates result from the matrix multiplication a^6 c 

'VV/aT 1 



V 



<r 



-v 
-v 



-VV(a 

2|V| 2 7 
7 2 + £ 2 



More generally for any n we have 
for i 7^ j and 



( N o) kk = -J2^ kl 



(12) 
(13) 



This is just the network described in [3J and the introduction. 



3.4 Including re-emission, recombination and trapping 

The population decreasing effects of re-emission, recombination and trapping can all be described by the rates Kk 
of the diagonal anti-hermitian operator 

n 

k=i 

included in our general master equation ([2]) . The contribution to the rate of change p is easily calculated 

-i{A,p} = -^2^{K k + > 



and M becomes 



M 



C\ —a 1 
a b + c 2 
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N = a f (6o + c 2 ) 1 a + ci 
N = a f (6 + c 2 )~ 1 a + c\ 
N k = a+(6 + C2)- 1 (-v(b + C2)- 1 )" a 



(14) 



which also hold with the replacements a — > a, b — > b and v — > while leaving c\ and C2 unchanged. 
The rates in N can again be calculated directly 




for k 7^ I and 



( N o)kk = ~Kk -J2^ kl ■ 



l=£k 



3.5 Numerical simulations 

According to the last two subsections, network Nq is easy to calculate directly, while network N and any fc-site 
contribution Nk can be formed from the general definition of a, bo and v (see Appendix (jEJ)) which can be somewhat 
tedious. However, there is another approach related to numerical simulations. When running numerical calculations 
to simulate a complex master equation ([2]) on a software like Octave or Matlab, the need to convert the equation 
to the form 



with a real M arises in any case. This can be done as we describe it in \2.2\ or more easily - because we have the 
help of a computer - by defining an orthonormal density space basis. For example set 



for k < I. Then matrix M can be formed by applying the master equation to those vectors and finding their 
coordinates. The following gives the population space block of M 



where M is the superoperator formed by the RHS of the master equation ([2]) . Once the entire real matrix is found 
it is cut into population and coherence blocks 



Hence, if one has already calculated M in order to simulate a quantum network, it only takes a few steps to find 
the kinetic network approximation N. 

4 Preliminaries 

In this section we give some definitions and conditions. The conditions allow us to infer basic facts about the spectra 
of the operators N , N and M, which are required for all our bounds in Sections ((5|), © and (JSJ) . 



P = Mp 



\k)(k\ 

(\k)(l\ + \l){k\)/V2 
(-i\k)(l\+i\l)(k\)/V2 




\mcp m C c J 

and a generalized kinetic network of the same form as N is calculated as 

N = mpcmccfncp + mpp . 



M = ( mpp mpc ] 
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4.1 Norm 



Because for our bounds of the relaxation time we remove all population decreasing effects all evolutions M, No, 
and N leave the total population invariant. Therefore we split up the space of populations P. Set 

e = (l,l,...,l)t/neP 

the equal population vector. As we will prove in Proposition 01 as long as the network meets certain conditions, 
both quantum and kinetic evolutions will tend to e for any initial condition with total population 1 . Consequently, 
we are only interested in the properties of our matrices in the space of population inequalities 



I = e T = (v 



This is reflected in the norm we use, defines as follows. For A : X\ — > X2 where Xi and X2 are equal to / or C 
we define the operator norm as 

11 411 !Ld«jla 
\\A\\ = sup 

vex 1 IMI2 

where |J-|| 2 is the Euclidean norm. Hence, from now on, we think of our matrix blocks as 

a:I ->C 

Note that ||a|| is the same if we maximize over I or P because ae = 0, and that ||a|| = || || . Also, according to 
Proposition [4] Ao < on / and therefore A^ 1 is well-defined. The same holds for N. Define 

M=ii^ _i ir 1 

to be the eigenvalue closest to in A^ in /, define /zq the same way for Nq. 



4.2 Conditions 

For all our following bounds we have a set of conditions. 

• First, we require that the network is connected, in the sense that any two sites k and I are coupled, at least 
via some intermediates, i.e. for some integer p > there are sites rrij, j = 1 .. .p such that the product 

^mi Knim2 ■ • • * nip— imp 'hp j 

is nonzero. This condition ensures that all sites can exchange population and the evolution ultimately con- 
verges to e. 

• Second, we require all the site dephasing rates to be strictly positive, 7^ > 0. This condition is essential for 
our approximation, as the coherences need to decay for the evolution M to become non-oscillatory. Notice 
that the limit 0r -1 — > does not require that the dephasing rates get larger, but they will be much larger 
than the magnitude of eigenvalues of or Nq, the population decay rates, because T 3> Q 2 T~ 1 . 

• Finally, we require that the Vjy are real. This ensures that A^ is symmetric and has a real spectrum (see 
Proposition [3]) , which allows simpler bounds in our proofs. While No is always symmetric, we first compare 
the evolutions of M and A", and then the evolutions of A" and Nq. Therefore we require this condition for 
both Nq and N. We are confident that our methods would extend to the case of complex but for the sake 
of clarity we restrict ourselves to the simpler case. 
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4.3 Inverse bounds 

Or proofs consist mainly of using the following two bounds on the inverse on different parts of resolvents. 
First, consider the Taylor series of the inverse close to 1, which for real numbers x gives 

|(1 + a;)" 1 - 1| <2|x| 

for |x| < 1/2. This is readily translated to a bound for operators 

\\{A + B)- 1 - A- x \\ <2\\A- 1 f\\B\\ (15) 

for ||B|| < | H^ir 1 . 

Second, if A < —c < is a negative definite, self-adjoint, finite dimensional operator and z G C with Rez > 
then 

IKz-a)- 1 !! < c- 1 (i6) 

and 

Wiz-Ay^Klzr 1 . (17) 

these two bounds follow from the fact 

\\(z - A)~ 1 \\ = max {|A||A G Spec (z - A)^ 1 } 
= max { I (z — A) 1 1 1 A G Spec A} . 

4.4 Spectral properties 

The following Proposition gives some basic facts about the spectra of the kinetic networks N and Nq. We will use 
these properties for the proofs of our bounds. 

Proposition 4. The matrices Nq and N as defined in \S.3\ and \3.1\ have the following properties 

1. Nq is real and symmetric. 

2. N is real. 

3. If the interactions Vu are real then N is symmetric. 

4. N Q e = Ne = 

5. If > and the network is connected network then Nq < on I . 

6. For 0r _1 small enough, Nq < — and N < — /io/2 on I. 

Proof. 1. These properties follow directly from the form in (fT2|) and (fT3| . 

2. N is real because it is a product of a, & _1 and which are also real. 

3. If Vm is real then a (see (jTUJ) ) is real, so 

N = fi T 6" 1 a 

Furthermore, b T = b (see Appendix [A]) therefore 

N T = ~a T (b- l ) T a 

= ~a T (b T )~ X a 
= N . 

4. From (|10[) it is not hard to understand how a looks for any n. One sees that the two rows for the coherence 
between sites k and I have exactly two non-zero entries, the first has Vki and — Vm, and the second has Vjy and 
—Vki- Therefore ae = and so Noe= Ne — 0. 

5. For v G / we have J2 k Vf. = 0. Now 

v^N v = ~y] fi k i(vk - vi) 2 

k<l 

< o 
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The condition for equality is as follows. Because 7^ > we have 



V k i ^ o m ± . 

Hence, because the network is connected, we have — vi for all k and I and with J^ fe Vk — it follows that t? = 0, 
thus No < on J. 

6. Because /ito = HA^ 1 )! an d Ao < we have Ao < — Mo on J. Note that /i cx 2 r _1 can grow as 0r _1 gets 
small, so [i — fiQ can grow in absolute value, however, as we now show the spectra of iV and Nq approach each other 
relative to their "size" 

H - fi < Mo 

We bound the distance of TV and Ao with the inverse bound. For 0r _1 small enough we have \\v\\ < | H&o" 1 !! 
and we can apply (|15p on 

||(6 + ^)- 1 -& - 1 || <2||6 - 1 || 2 || i ,|| . 



Now 



||JV — JVbll = ||a f (b- 1 -^ 1 ) 



= IH| 2 2||6 - 1 || 2 |H . 

So, the distance of N and A^o is proportional to O 2 r _2 0, and the eigenvalues in No and N - in particular /j,q and 
/i - are proportional to 2 r _1 . Comparing the two gives 

e 2 r~ 2 e < e 2 r~ 1 

because 0r _1 <C 1. That means the eigenvalues are approaching each other relative to their magnitude, in particular 
N becomes negative definite like No, and 

Im-moI 



Mo 



0. 



Now it immediately follows that 



N < -fi< (i /2 
N < -Mo < M/2 ■ 



□ 



5 Bounding relaxation time error 

We now give an explicit definition of relaxation time and the norms we use to control it. Then we derive bounds 
first comparing the quantum network M to the kinetic network N, and then comparing the kinetic networks N and 
iVo- 

As a simple check of sanity consider the following. If we scale r cx s and cx s then also M, N cx s and time 
scales inversely At, t cx s _1 . Therefore the relative error Ar rc i = Ar/r stays unchanged and we expect bounds 
in terms of positive powers of 0r _1 . Our two bounds show exactly this behavior. The approximation of A^ to 
M is proportional to 2 r -2 , while the approximation of No to A^ is proportional to 0r _1 , combining the two 
approximations it follows that the approximation of Nq to M is also proportional 0r _1 . 

Note that all the results in this Section require the conditions in 14.21 

5.1 Relaxation time 

By Proposition (J3J), the eigenvalues of N and Nq on / are all negative for ©I 1 " 1 small enough, so for any initial 
distribution po G I 

e Nt po 
e Not p 
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for large t. We can integrate 



e m dt = N~ 1 



and applying the operator norm maximizes the relaxation time for the kinetic network N over all possible population 
inequalities po £ I, set 

t = /i- 1 = IliV- 1 !! = II (aV^)- 1 ]! 



dt 



and in the same way we define tq = fj, 1 for the network No. 

We define the error in relaxation time as the relaxation time difference maximized over / 



At = 



Te Mt T i 



r.Nt 



dt 



Hence, bounding At means controlling the worst possible error in relaxation time when approximating M by N. 
The relative error is 

AtYci = At )t 

notice that we compare the worst possible relaxation time error to the longest possible relaxation time, those two 
do not necessarily occur for the same initial condition. We define Atq and A To, re i in the same way, comparing N 
and Nq. 



5.2 Resolvent difference 

Converting the operator for the relaxation time error we get 

> 

Te Mt T j _ e Nt dt = TM -1 T 1 _ N -l 

1 



— (f)- 

2wi j z 
1 / 1 



T- 



z-M 



z-N 



dz 



2ni J z \z — a< (b — z) 1 a z — cvb 1 a 



dz . 



where the complex integration follows a contour surrounding both SpecM and SpeciV. Define S(z) to be the 
difference of the two resolvents 

S(z) 1 



We now seek a bound on 



z — a'(b — z) 1 a 



Te Mt T^ - e Nt dt 



2m T ■ 



(18) 



5.3 Comparing the relaxation time of M and TV 

When bounding second order terms with the inverse bound we encounter 

« = llall 2 !!^ 1 !! 2 

and kq for the corresponding terms with bo instead of b. Notice the scaling behavior /x,/io oc 2 T -1 and k, ko oc 

e 2 r- 2 . 

We will change the contour integration in f| 18 j) to be along the imaginary axis z = iy for y G R. We prove the 
somewhat technical bounds on S(iy) in Lemma [TT1 in Section [S] 



Theorem 5. If QT 1 is small enough then 



At = 



Te Mt T^ - e m dt 



7T 
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where j3 > is the scaling independent constant from Lemma This gives a bound on the relative error 



Ar re i = At/t < -k(1 + 0) = k 2 e 2 T- 

IT 



where k 2 is scaling invariant. 



Proof. We set the integration contour in f| 18[) to be along the complex axis z = iy for y € R with y going from —R 
to +R. We close the contour to the left in the half plane of negative real parts along a circle of radius R. According 
to Lemma 1111 S(z) has no poles with Re z > and so all poles lie within this contour for R large enough and 
Or -1 small enough. As R tends to infinity the integrand behaves like so the half-circle does not contribute to 
the integral. We can therefore change to complex integral to an integral in y over all of R 



z m dt 



r~ / —S{iy)dy 
2tt J r iy 



Now split up the integral into two regions \y\ < fi and \y\ > fj, and then use the corresponding bounds from 
Lemma [TT1 Choose 0r _1 small enough so that /i < a and use part 1 of the Lemma to bound 



/ —S(iy)dy < / — ■ An^" 2 \y\dy 
J-n W J -a \y\ 



and use part 2 of the Lemma to bound 



-S{iy) dy 



ii \y\ 



i 



< / — • A[3k\ V \- x dy 

Jfi \y\ 



Adding the two bounds gives the result 



< 



/ — S{iy) dy 

2m J R iy 



< —Sk/*- 1 (1 + 0). 



□ 



5.4 Comparing the relaxation time of N and iVo 
Theorem 6. // 8r _1 is small enough then 



An = 



,Nt „N t 



dt 



<4^ 2 ||^|| 



where fi and k can also be replaced by ^ an d k . This gives a bound on the relative error 

An ir ei = Ati/t < Ak^i' 1 \\v\\ = ^Gr^ 1 

where k' 2 is scaling invariant. 

Proof. In this case we don't need to bound the resolvent, instead we can evaluate the integral 

e Nt -e Not dt = N- 1 -N^ 1 . 



We use the inverse bound (|15[) twice. First, because < \ ||b 1 || 1 as long as QT 1 is small enough, we can 
apply the bound on 



I (6 i^) 1 & 1 II <2||6- 1 |r||r/|i 



(19) 
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Now, apply the bound again with A = N and B = No — N . The condition for B is 

\\B\\<\\af\\(b-v)-i-b- x \\ 

= 2 «IMI 

< 2 = 2n 

where we used ([!§]) in the second step. The last inequality is again achieved for 8r _1 small enough because the 
two sides scale like 

e 2 r~ 2 e < e^- 1 . 

Now it follows that 

ii 1 - TV- 1 ! ^ 2 ii^ii 2 = 4k ^ 2 hi 

as claimed. By switching the role of b and bo we receive the corresponding bound with kq and /iq. □ 

As a corollary we receive a bound on the relaxation time difference between the fully quantum mechanical 
evolution of M and the simple kinetic network evolution of No. 

Corollary 7. If Or -1 is small enough then for some scaling independent constant ki 

Aro.rei < her- 1 

Proof. According to Proposition [4] we have — //o| /a*o ~~ ► 0> an d therefore there is a c such that 

c > t/tq 

for Or -1 small enough. Then with Theorems [5] and [5] we have 

AT , rol = Atq/tq 

< (Ar + Ari)/r 

< c(At + An)/r 

< c(Ar re i + An )re i) 

< fcier- 1 . 

□ 

6 Bounding evolution error 

In this chapter we bound the difference of time evolution operators for M, N and A*o. Our error bounds looks as 
follows 

|| e Aft_ e ATt|| < e -pt/2. X 

where X is proportional to Q 2 T~ 2 up to a logarithmic term, and proportional to Or -1 if N is replaced with N . 
The logarithmic term appears due to intermediate times. It seems the integral over time performed in the last 
chapter seems to have conveniently guided us around that logarithm. As for the time dependence, using a shifting 
integration contour might give a bound like e~ M */rf, but a better control of the spectrum would be necessary to 
shift the contour close to — fi for long times. 

As in the last chapter in l5.2[ we write the evolution difference as a complex integral before we prove bounds 



\\ Te Mt T ] _ e m\\ J_I e zt s{z)dt 
2m J 

Note that all the results in this Section again require the conditions in 



(20) 
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6.1 Comparing the evolution of M and iV 

We will change the contour integration in ([20]) to be parallel to the imaginary axis z = iy — (i/2 for y £ R. With 
this choice the exponential in the integral yields exponential decay at rate /k/2. Again we give the technical bounds 
on S(iy — n/2) in Lemma 1X21 in Section |H1 

Theorem 8. If QT^ 1 is small enough then for all t > we have 

\\Te Mt T* -e m \\< e - "* • fc 4 e 2 r~ 2 (i + k 5 in e _1 r) 

where fc 4 and k§ are a scaling independent constants. 

Proof. We set the integration contour in (|20[) to parallel to the complex axis z = iy — fx/2 for y G R with y going 
from —R to +R. We close the contour to the left in the half plane of negative real parts along a circle or radius R. 
According to Lemmas [Til and [i"2l S(z) is bounded for Rez > —(J./2 and hence has no poles. Therefore all the poles 
lie within the contour for R large enough and 0r _1 small enough. As R tends to infinity the integrand behaves 
like p-e Rcz so the half-circle does not contribute to the integral. We can therefore change to complex integral to 
an integral in y over all of R 



|Te Mt Tt - e 



Nt I 



/ eto-^Sfa-vWdy 



Now split up the integral into three regions with \y\ in the intervals [0, fj], \pt, a] and [a, +oo) and then use the 
bounds from Lemma [121 Choose 9r _1 small enough so that /x < a and use part 1 of the Lemma to bound 



'to-f/WSiiy- (x/2)dy 



< e-^/ 2 16«: 



nfi/2 

/ l6K(j,~ 2 \iy - n/2\ dy 
Jo 



and 



n/2 



e (fv-M/a)*s(iy - n/2) dy 



not, 

/ An-\y\- 2 \iy- n/2\dy 
/ An-2y~ 1 dy 

Ju/2 



'A»/2 

e-^/ 2 8Kln(2a//i) 



and we use part 2 of the Lemma to bound 



< e-^l 2 



4|y|- 2 ||a|| 2 (6 min - Ai /2)- 1 ^ 



<e-^ 2 Aa- l \\a\\ 2 {b ljAn -^/2)- 1 . 

Adding the three bounds gives the result 

\\Te m Tl - e Nt \\ < e"^ t/2 4 (4k + 2 K ln(2a/^) + a" 1 ||a|| 2 (6 min - . 

The middle term of the parenthesis has the worst scaling behavior 

2k ln(2d//i) cx 9 2 r~ 2 In 9 _1 r 

while the other two terms scale like 9 2 T~ 2 . Therefore there are some scaling independent constants fc 4 and k$ such 
that 

\\Te Mt T ] - e Nt \\ < e^ t/2 • k^T' 2 (l + k 5 In _1 r) . 

□ 
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6.2 Comparing the evolution of N and iVo 
Theorem 9. If QT^ 1 is small enough then for all t > we have 

|| e "t_ e JVot|| < e -M*/2 . ^er- 1 

where is scaling independent, and where \i and k can also be replaced by [1q and Kq. 
Proof. We are bounding the integral 

^-ie zt S{z)dz 
2m f 

with resolvent difference 

S(z) = — — . 

K ' z-N z-N 

We use the same contour as in Proposition [5J z — iy — fi/2. According to Proposition 01 all poles of S(z) lie within 
this contour when Or -1 is small enough and R is large enough. Because of the e zt factor and T(z) tending to zero, 
the integral over the half-circle tends to as R becomes large. 

We bound S(z) in much the same way that we bounded S(z) in Lemma 1121 however, the procedure is more 
straightforward. Set 

X = {b-is)- 1 -b' 1 , 

because |M| < h ll^ 1 !' 



we can use the inverse bound (|15[) 

.-111 2 : 



\X\\ < 2 W 



HI • 



Now rewrite 



S(z) = (z- a^b^ay 1 + (z - atb^a - a^Xa)' 1 
For any z with Re z = — (a/2 and for 0r _1 small enough we have 



\atXa\\ < 2 ||a| 

1 M 



Ml 



< — \\z — a)b 1 a|| 
and so we can apply (fT5|) again 

S(z) < 2||(z-a t 6a)- 1 || 2 \\a*Xa\ 
< 4\\(z - a^ba)' 1 ^ k\\p\\ . 
Now we apply inverse bounds and (|T7)) to receive the bounds 

S(z) < 16/i~ 2 K||^|| 



S(z) 



< 4 \z \ k M 



as long as Rez = — /i/2. 
To estimate the integral 



1 

2^ 



e^-^Tiiy - fi/2)dy 



we split it into the two regions \y\ < n/2 and \y\ > /i/2. Use (|2"Tj) to bound 

rf i/2 



n/2 



< e-" t/2 -fx- 16^ 2 k|M| 



and use (1221) to bound 



At/2 



,(iy-v./2)t§( iy _ n/2)dy 



< 2e-^ t/2 ■ Sn^nl 



(21) 
(22) 
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Adding the two bounds gives 

|| e M_ e w t|| < 32e-^ 2 f i- 1 K\\u\\ 
< e~^ /2 ■ keOT- 1 . 

where ke is scaling independent. The whole proof works just as well when exchanging (i with fj,Q, k with kq giving 
a similar bound. □ 

As a corollary we receive a bound on the decay time difference between the fully quantum mechanical evolution 
of M and the simple kinetic network evolution of Nq. 

Corollary 10. If 0I 1-1 is small enough then for all t > we have 

||Te Mt rt - e w °*|| < e~^ 2 ■ hOT- 1 
where k% is a scaling independent constant. 

Proof. The bound follows from Theorems [5] and |H] and the fact that 

e 2 r 2 ine _1 r < er^ 1 

for er- 1 < 1. □ 
7 Applications 

The rate of direct population exchange 



Hkl 



2\V k i\ 2 Jki 
lit + El 



determines the strength of the link between sites k and I for the network Nq. Because of our condition that 7^ > 0, 
the network topology is fully determined by the Vki, but the strength of the links is also affected by jk and E}.. 

As applications, we consider two idealized networks. The first is a highly connected network where all sites are 
linked, the second is a circular chain where where only nearest neighbors are linked. We numerically calculate the 
relaxation times for the networks M , Nq and N and compare the relative errors. Then we compare these networks 
to randomized networks with the same network topology. We also discuss the dimension dependence of our bounds 
from Sections [5] and again compare it to numerical simulations. All the simulations agree with our bounds, but 
they show much room for improvement when considering large dimensions. 

Finally, we discuss the FMO-complex and our model for which some results were already shown in the intro- 
duction in Figure [TJ 

For clarity of notation we recall that At, Atq and Art are relaxation time differences between the network pairs 
M — N, M — No and N — Nq respectively. This only makes the discussion more precise, while generally Atq and 
At} show the same dimension and scaling behavior, with small corrections to constants. 

7.1 Highly connected network 

Consider a highly connected network 

v k i = e 
E k -0 
7fe = r . 

In Figure[2]we made a plot for the computed relative relaxation time differences Ar re i and Aro, re i for different Or -1 
with the initial state localized at site 1. Both axes plot logarithms, hence a straight line with slope n represents a 
(Gr -1 )' 1 proportionality 

The difference Ar rc i is too small to show any clear behavior. The difference Aro. re i is linear with slope approx- 
imately 2, hence the approximation is better than the slope 1 expected from Theorem [6] In the same figure we 
compare our idealized network to random networks where all Vki are chosen randomly between and 9 and all Ek 
are chosen randomly between and T, hence they have the same topology. The magnitudes of the errors are similar 
for the range considered, but the slopes are different. All the samples show an error slope of 1 for Atq^oI, while 
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-5 -4 -3-2-10 1 2 

Figure 2: Relative error for the highly connected network 



the error slope for A-r re i is varying, but in most parts steeper than the slope of Ato, rel- This behavior is closer to 
the behavior expected from our bounds. Generally, the agreement is about six orders of magnitude better for the 
network N than the network No- 

For the ideal highly connected network we derive the quantities used in Theorem [5] and |5] analytically in Ap- 
pendix The resulting bounds are 

Ar rc i < Cl n0 2 r- 2 
Ar^rei < c 2 rt6r~ 1 

for dimension and scaling independent constants c\ and c 2 . The simulation of M has a relatively high error and 
becomes slow very fast as n gets larger. Hence, we can only get meaningful results for An !re i, the relaxation time 
difference of networks N and No. The result in Figure [3] actually shows that the difference increases with slope 
2 or proportional to n 2 . The reason is that in Theorem [5] we have the condition||^|| < ^ where the LHS 

is proportional to n and the RHS is constant (also discussed in the Appendix). If we increase the dimension at 
constant scaling, this condition and our bound break down. To still get a bound for large n we would need to 
readjust the scaling. 



7.2 Linear network 

Assume the sites are positioned on a circle and only nearest neighbors interact with strength 

V kl = 

where we use the equivalence n = 0. Further "fk = T and such that E^i = TE when \k — 1\ = 1 which is possible 
for n even. 

In Figure |4] we made a plot of the computed relative relaxation time differences Ar rc i and Ato, rci for different 
0r _1 with the initial state localized at site 1. Interestingly the quality of approximation by Ao is improved over the 
highly connected model, while the quality of approximation by N has decreased. Also, both models show the same 
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Figure 3: Relative errors between N and Nq for the highly connected network and the cyclical chain with increasing 
dimension and = 0.01 and T = 1. 

slope of about 2. We compare the ideal chain to random chains for which the 14/ that equal in the idealized case 
are instead chosen randomly between and 0, and all Ek are chosen randomly between and T. We get essentially 
the same behavior with all slopes being 2. That hints at a possible improvement of our bound in Theorem [5] in the 
case where the network is a chain, improving the proportionality from Or -1 to 2 r~ 2 . Generally, the agreement 
is about five orders of magnitude better for the network N than the network Nq . 

As in the last section, we can derive the necessary quantities for our bounds and get 

Ar rcl < c 3 2 r- 2 7i 2 
ATi, rc i < c 4 0r~ 1 n 2 

for dimension and scaling independent constants C3 and C4. This time the conditional j < i 1 does not 

break down and the bounds hold for large dimensions as well. The n 2 terms are due to the lowest eigenvalue of 
Nq being proportional to iiT 2 . This is a weakness of our strategy to use the operator norm for our bounds. Better 
bounds should be possible when only considering localized exciton as initial state. This initial state would the a 
superposition of all the eigenstates on No, and the average relaxation time would enter the bounds, instead of the 
longest relaxation time (the smallest eigenvalue of Nq). 

As above we skip the simulation of M because the error is too large, and consider Ari, re i only. The result 
in Figure [3] shows that the difference seems to approximate a constant value for larger dimensions. So, both our 
bounds could be improved for large dimensions. 



7.3 The FMO-complex 

The FMO-complex is pigment-protein with trimer structure. Each monomer contains seven bacteriochlorophyll a 
pigments that capture and transport light. The excitons start out at site 1 or 6 and the trapping occurs at site 3 
[1], we set the initial state to be 

#, = (1/2,0,0,0,0,1/2,0)*. 
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io gl0 i/r -1 

Figure 4: Relative error for the circular chain 



We use the same numerical values as with interactions and energies from [¥]. The system Hamiltonian is 



H + V = 
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with all the numbers in cm 1 (or 2.9978 • 10 10 s 1 ). Exciton recombination at rate k — Ins 1 and reaction center 
trapping at rate k 3 = lps -1 enter the anti-hermitian operator 

A = - t -^K\k}{k\ + K3 \3}{3\\ . 

We use the same dephasing rate for every site 7^ = 7, and vary 7 from 10~ 3 to 10 5 cm _1 . Efficiency is calculated as 

/>oo 

/ = k 3 / p 33 (t) dt 
Jo 

we calculated / for the three models in Figure Q] Peak efficiency is reached for 7 ss 170cm -1 close to the average 
energy gap along the chain which is 146cm -1 . The approximation N has less than 1% error, even for the lowest 7 
used, and the approximation Nq gets below 1% error for 7 « 2cm~ 1 . Comparing this to our bounds we have 

||a|| = Hatjl = 215cm -1 

and for large 7 
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The numerical factor /3 is changing because of the changing ratio between energies and dephasing, for large 7 
however it is approximately equal to 100. Hence, our bound becomes 



Ar rcl < 100 (215c 



The 1% error margin is reached only when 7 = 21500cm -1 , so our numerical factors could certainly be much 
improved. But this is not unexpected, since our main goal was to find the leading behavior in Or -1 . 
We give Nq for maximal transfer efficiency 



jV (7 = 170cm- 1 ) 



It is interesting that the rate between sites 2 and 3 is actually smaller than the rate between sites 2 and 6 even 
though I V23 1 > I V26 1 ■ The reason is the large energy gap between sites 2 and 3 of 420cm _1 while sites 2 and 6 have 
an energy gap of 60cm _1 . However, the values for site energies are still up to some debate [JJE], and small changes 
can easily turn this behavior to the opposite again. 
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8 Resolvent difference bounds 

The following three Lemmas are the main technical parts of our bounds. They all consist of bounding the operator 
norm of the resolvent difference 

S( Z ) = 1 

z — a'(b— z) 1 a z—a'b 1 a 

for different values of z. Conceptually the bounding procedure is simple, we only employ the inverse bounds 
introduced in 14.31 Loosely speaking, if \z\ < T we can expand (b — z)~ l and then the two terms in S(z) only 
have a small difference in the denominator, so, using another inverse bound, they almost cancel. If \z\ > T then 
\z\ ^> 1 1 ci+ £>~ 1 a 1 1 an d we can directly use the second step from the case \z\ < L. 

Of course we also have to keep in mind where the poles of S(z) are. According to Proposition [H (z — A) -1 has 
poles on the real axis below —fi which move according to the scaling 8 2 L _1 . On the other hand (z — a^(b — z)~ 1 a)~ 1 
has poles close to the poles of (z — A) -1 that approximately cancel each other, but it also has poles close to the 
eigenvalues of b which are approximately etij = —Jij +iEij and ciy, scaling like T. Comparing the two sets of poles, 
the 6-poles are much further to the left (negative real values) than the TV-poles because V 3> 2 r _1 . Our lemma 
steer clear of this poles by keeping Rez > — /z/2. 

Lemma [IT] contains bounds for Re z > which on the one hand ensures there are no poles on the right side of 
the complex plane, and on the other hand we use the bounds for z = iy to bound the relaxation time. Lemma [T2l 
contains bounds for the region — ji/2 < Rez < the bounds are derived in a similar fashion as in Lemma 111 [ but 
there are some additional complications. 

8.1 Bounds in the right half plane 

Lemma 11. If QT^ 1 is small enough and Rez > then S(z) is bounded by 
1. \\S(z)\\ < 4K/i _2 |z| if \z\ < a, where a oc T depends on a and b, 

2- HiS^z)!! < 4 / 3k|z|~ 1 for any z with Rez > 0, where (3 is a scaling independent constant depending on a and b. 
Proof. 1. Assume Rez > and |z| < a oc T, where 

a = min Her 1 ! 1 ,^n^ 1 fi\. (23) 



Set 



X = (b - z)- 1 - b- 1 
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and because \z\ < 5 II 6 1 || 1 we can use (fT5j) and have 

\\X\\ <2||&- 1 f|z|. 

Rewrite 

S{z) = {z- cjb^a - a^Xa)- 1 - (z - aV^)" 1 
To use (fl"5|) on this expression notice that 

, 1 -1 

\z\ < —K fl 



and therefore 



i^Xa\\ < 2k\z\ 



1 

< -Lb 

~ 2 P 

< - 1|( 2 -aV^r^r 1 

where P^|) was applied in the last step, using the fact that a%~ 1 a is self-adjoint from Proposition @] This is just 
the condition for the bound 

\\S{z)\\ <2\\(z- aV 1 ^)- 1 !! 2 W^XaW 
< 4k||(z- a + 6 -1 o) -1 |[ 2 |z| 
again using (fT6]l and also (fl~7f we get the bounds 

\\S(z)\\ <^- 2 \z\ 

\\S(z)\\<4^z\- 1 (24) 
for \z\ < a. The first bound is bound 1 of the Lemma, the second bound will be used below. 



2. We now derive a bound when \z\ > a and Re z > 0, we will combine it with (|24[) to receive bound 2 for all 
If 9r _1 is small enough then we have 

Hawaii < - < -|z| 
II 2 2 

||a f (6- z) _1 a|| < - < -\z\ . 
H v , 2 2 

Where the latter inequality uses the fact that the spectrum of b approaches the spectrum of bo as 0L~ X becomes 
small, and the spectrum of bo, which is —7^ ± iE%j, has negative real part — 7^ < 0. The last two inequalities are 
the conditions to use (|15|) and get the two bounds 

|| (z - a\b - z)- 1 ^- 1 - z- 1 !) < 2|z|- 2 ||o + (6 - 
||(z - ot^a) -1 - z" 1 !! < 2|z|~ 2 lla+b-^ll 

set 

6min = min {| Re A|| A e Spec 6} cx T (25) 
the closest any eigenvalue of b gets to the imaginary axis. Then ||& _1 || < ^min an d \\(b — — ^min so 

||5(z)||<4|z|- 2 || a || 2 d- 1 . 



Comparing to (|24l) with 
we have 

for \z\ > a and therefore 



(3 = max 1 1 , 1/ (abnxin \\b 1 1 1 | cx 1 
4|«|- a ||o|| 2 6^ < ;SF - ^j^l- 1 



<^k\z\~ x (26) 
for all z with Re z > 0. Which is bound 2 of the Lemma. □ 
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8.2 Bounds parallel to the imaginary axis 

The following Lemma establishes bounds along the imaginary axis z = iy — p. These bounds are used to prove the 
evolution bounds. 

Lemma 12. // we choose QT^ 1 small enough then for < p < fj,/2 the resolvent difference S(iy — p) is bounded 
by 

1. \\S(iy — p)\\ < \Qnii~ 2 \iy — p\ and \\S(iy — p)\\ < AK\y\~ 2 \iy — p\ if \y\ < a, where a cx T depends on a and b, 

2. \\S(y)\\ < 4\y\- 2 \\a\\ 2 (b min - ft)- 1 for \y\ > a with b min cx T. 

Proof. We proceed almost identically as in the proof of Lemma [11] using the inverse bounds [15j (fl"6)) and (fl7)) for 
the same parts of the resolvent terms. 

1. We use the a from (l23l) to define 



1 ||. — x II — i 1 -i 
a = mm 1 1 II > g K M r" ~ M 

notice that the scaling a cx T is only approximate and that 8r _1 needs to be small enough such that a > 0. Now 
require \y\ < a cx T. Set 

X = (b-iy + fi)- 1 -b- 1 



and because we have 



\iy-P>\ < 1 2/1 + A 



we can use (fT5[) to get the bound 
Rewrite 



< J II 6 1 1 + A 

<-lP-T' 



\X\\ <2\\b- 1 \\ 2 \iy-fi 



S(iy — p) — (iy — p — a b 1 a — a) Xa) 1 — (iy — p — a)b 1 a) 1 
To use (fT5j) on this expression notice that we have 

1*2/ — Al <\y\ + P 

< [ -k~V — M ) + A 



1 

8 ^ 

and therefore 

1 1 a < 2n\iy — fi\ 
1 

< -a 
~ 4^ 

< o ||(*2/ — M — o t 6 -1 a) _1 || 

where (JTHJ) was applied in the last step, using the fact that p + a^b^a < — /i/2 from Proposition @] This is just the 
condition for the bound 

\\S(iy - p)\\ < 2 \\(iy - p - alb^a)- 1 ]] 2 \\atXa\\ 
< 4k II (iy — p — a^6 _1 a) _1 II \iy — p\ 



25 



again using (fTtl)) and also (fTT)) we get the bounds 

\\S(iy-fi)\\ < 16Kn~ 2 \iy - p.\ 
\\S(iy-p)\\ <4K\y\- 2 \iy-fi} 

for \y\ < <x These are the bounds in part 1 of our Lemma. 

2. We now derive a bound when \y\ > a. If 0r _1 is small enough then 



a f (6 - iy + p)~ l a + p\\ < ^ < ^|y| . 



The last two inequalities are the conditions to use (fTS"]) and get the two bounds 

|[ (iy - /* - a f (6 - iy + m) -1 ^ 1 - (*» - < 2|yT 2 ||ot(6 - iy + /i)" 1 

|| (iy - /2 - aV 1 ^- 1 - (iy - /i)" 1 )) < 2|y|~ 2 ||a t 6 _1 a|| . 

Use & m in from PS]), giving 



a 



{b~iy + p) 
(6 - iy + A) 




A) 
A) 



-l 



-l 



and so 



\\S(iy-fl)\\ <4|y|- 2 ||a|| 2 (6, 



A) 



-1 



for |y| > d, which is the bound in part 2 of our Lemma. 



□ 



9 Conclusion 

We studied to kinetic networks that approximate the energy transfer in a quantum network subject to dephasing. 
The first network N derives its rates only from nearest neighbor interactions, while the second N includes higher 
order corrections. We proved that the relaxation times are proportional to 6L _1 and 2 L~ 2 respectively. Hence, 
the approximations are good if the interaction gets weak, or the dephasing and/or energy gaps get large. In the case 
of the FMO complex, both kinetic networks are good approximations in the regime of dephasing-assisted energy 
transfer. With simulations we found that the more complex kinetic network N provides approximations with a 
percentage error 5-6 magnitudes smaller than the simple kinetic network. 

The study of these approximations could be extended in several ways. First, one could study the higher order 
corrections involved in N. Second, when the interactions Vm are complex, N can be non-symmetric, meaning 
population exchange between sites is directed, this might relate to coherent cancellations along loops as mentioned 
in [3]. And finally, it would be interesting how our method of splitting population and coherence space to achieve 
kinetic network approximations could be generalized to other quantum networks and how it relates to existing 
models to approximate coherent evolution with incoherent statistical evolution. 
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A Three sites 

In the following we write out parts of the master equation §3§ for the case n = 3 and then derive the form of the 
matrix M. Then we explain how to generalize that form to higher n. For simplicity of notation we omit the scaling 
factors r and 0, until we reach a block matrix expression. First note that with a standard calculation one finds 
£(p) to decrease the coherences in the manner 



(£(p))ki = -IkiPki 
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where k ^ I and = 5(7^ + 7/) and (C(p)) kk = 0. This gives a diagonal contribution —Jm in the diagonal of the 
two rows corresponding to the real and imaginary part of ■ 
Now, we evaluate the commutator 



Ei 


V12 


Vu\ 1 pn 


P12 


Pl3 


V21 


E2 


V23 , P21 


P22 


P23 


Vsi 


V 32 


E 3 J \p31 


P32 


P33 



From the lxl entry we get 

Pll = -i(ElPu + ^12^21 + V13P31 - Eipu - V"2lPl2 - V31P13) 

= -i(Vi2Pu + V13P13 - V^i 2 pi2 - V^i3Pi 3 ) 
= 2Im(Vi 2/ oi2 + Vi 3 pi 3 ) 
= H-v^p\ 2 + n 2 p r i2 ~ W3P13 + Vl 3 pl 3 ) 
where superscripts r and i are shortcuts for real and imaginary parts, and from the 1x2 entry we get 

P12 = ~i(EiPi2 + V12P22 + Vi 3 p 32 - V 12 pn - E2P12 - V32P13) - 712P12 
= -i((Ei - E 2 )pi2 + V12P22 - V 12 pu + Vi 3 p 23 - V23P13) - 712P12 
with real and imaginary parts 

P12 = -VI2PU + VI2P22 - 712P12 + {Ei - E 2 )p\ 2 + Vi 3 p r 13 - V 2 r 3 p\ 3 + Vl 3 p r 23 - VI3P23 

Pvi = V{ 2 pU ~ V[ 2 P22 - (El - E 2 )p\ 2 ~ 712P12 + V&fis + y 23Pl3 ~ W3P23 ~ ^lSpIs ■ 

From these results we can read off lines 1, 4 and 5 of the following matrix and fill in the remaining lines in the same 
fashion 



M 






V2Vj 2 
-V2VI2 


-V2V{ 2 
V2V{ 2 


V2Vj 3 
-V2V? 3 


-V2V{ 3 
V2V{ 3 


V2Vj 3 
-V2V 2 \ 


-V2VJ3 

V2V 2 r 3 


-V2V? 2 


VWI2 




-712 


E12 


V 2 \ 


-V 2 \ 


Vh 


-v{ 3 


V2V{ 2 


-V2V{ 2 




— E%2 


-712 


V{ 3 


v 2 \ 


-V{ 3 




-V2V{ 3 




V2Vj 3 


v 23 


v 23 


-713 


E13 


v{ 2 


yr 
v 12 


V2V{ 3 




-V2V{ 3 


yr 
v 23 


v 23 


—Ei 3 


-713 


_yr 
v 12 


v 12 




-V2V 2 l 3 


s/2Vi 3 


v 13 


yr 
v 13 


-VI2 


yr 
v 12 


-723 


E 23 




V2V 2 r 3 


-V2V 2 r 3 


v 13 


V i 
'13 


-v{ 2 


v 12 


— E23 


-723 , 



V 

where we define Eij — Ei — Ea. Written as a block matrix 



M 



-at 
a b 



one can see the explicit form of the matrices a, and b. Remember that we also separated b into two parts. We set 
the 2x2-block diagonal that scales like T (the £!y and 7^ entries) to be bo and we set the block-off-diagonal that 
scales like (all the Vij entries) to be v. So b = bo + v. 

In 13.31 in ([H]) we defined a transformation U to diagonalize bo 1 if we extend this transformation to the entire 
space P C as 

U=l n ®U 

we can apply it to M directly and get 

/ 



M = U^MU 






-Vl 2 

V12 


-V12 
Via 


-Vi 3 
V13 


-v 13 

V13 


-V23 
V23 


-v 23 
V23 


V12 


-V12 




OL\2 




-iV 23 






-1V13 


V12 


-Vu 










iV 2 3 


iVu 




Via 




-V13 


-iV 23 




Q!13 




iVi 2 




V13 




-V13 




iV 23 








-iV X2 




v 23 


-v 23 




iVi 3 


IV12 




OL 23 






v 23 


-v 23 


-iVia 






-1V12 




OL23 
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row 


column 


entry 


row 


column 


entry 


kl 


km 


-iVim 


kl 


km 


iVi m 


Ik 


mk 


iVlm 


Ik 


mk 


-iVi m 


Ik 


km 


-iVi m 


Ik 


km 


iVi m 



Table 1: The non-zero entries of v 1 always I ^ m 
where ay = — 7y + iEy. This new matrix consists of the matrices a and bo also introduced in 13.31 




b = U^bU 
= b + v 

with v = WvU. The two kinetic networks are 

N = cJb^a 
N = tfb^a 

which also holds with all the tildes removed. 

It is straightforward to generalize the matrices a and bo to n > 3. Matrix a connects the population of site k 
to the coherences between site k and any other site I with strength Vki, and matrix bo is a diagonal matrix with 
entries ay and ay . A bit more complicated is the matrix v it is described in the next subsection. 



B General construction 

Here we give a description of how to find a, bo and v for general n. We number the n dimensions of population 
space P with k where k = 1, 2, . . . n and the (n 2 — n) dimensions of coherence space C with kl and kl where k < I 
are numbers from 1 to n. According to the order defined in 12.21 the first few dimensions of C are called 12, 12, 13, 
23, 23, 24, etc. . 

B.l Constructing a and b 

Matrix a is an n x (n 2 — n) complex matrix, with the only nonzero entries 

a-kM =Vu — —ak.ik 

hence in every column there are only (n — 1) nonzero entries. 
Matrix bo is diagonal with entries 

(b ) =~7M+iEki 

V / kl,kl 

(bo) = -Jkl - iEki ■ 

V / kl,kl 

B.2 Constructing v 

The matrix v = U*vU for any n is a somewhat complicated pattern of entries Vki, signs and complex conjugates. 
It connects coherences between sites k and I with coherences between sites k and m with the strength Vi m . Entries 
of v are only non-zero if one number of the two double indices match with further conditions on their conjugation. 
Table Q] shows the rules for the nonzero entries. 
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C Calculations for applications 
C.l Highly connected network 

Assume all sites are equally interacting, and have the same energies and dephasing rates 

v k i = e 

E k = 

ik = r . 

Then every column in a has 2{n — 1) non-zero entries all equal to 0. A simple calculation shows that 

a) a = 2nQ 2 (l n - nee') 
so for any v £ I we have a^av = 2Q 2 nv, hence 

||a|| = V2nQ. 

Obviously^ — —Tic an d ||^o 1 || = This gives k = 2n0 2 T~ 2 . Because 

a + &o l a = -r _1 a t a 

we have fio — 2n0 2 T _1 . Using u m uq we find 

• Pn^-iir 1 1 -i 

a = mm < — o , — k u 

1 2 H H '4 p 

. [1 12n0 2 r- 1 
= mm < — r, — = 

{2 '4 2n9 2 r- 2 

= r/4 

and so 

P = max jl, aT 1 
= 4 

and with Theorem [5] we get the bounds 

At < *r-i 
Ar rol < — n6 2 r- 2 . 

To get the bound on Ari )re i we also estimate \\v\\, we use the fact that each column and row of v has (n 
nonzero entries and so 

\\u\\ > vnQ 

for a scaling and dimension independent constant v. Then Theorem [5] gives the bound 

Ari, re i < AvnQT- 1 . 

The condition for this bound is 

\n<\\\b- l V 

the LHS is bounded from below by vnQ and the RHS is constant, so the condition does not hold for large n. 
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C.2 Circular chain 

Assume the sites are positioned on a circle and only nearest neighbors interact with strength 



Vi 



kl 



e |fc-z| = i 

else 



where we set equivalence n = 0. Further 7^ = Y and Eu such that Eui = YE when |fc — l\ = 1 which is possible for 
n even. 

Now, the column for site A; in a has only 4 entries, two each for the coherences with k — 1 and fc + 1. We calculate 

' 46 2 fc = I 

(aU) M = { -29 2 |fe-i|=l 

else 

So = 89 2 and ||a|| = \/8Q, in particular there is no n dependency. Also ||6q || = 1/Vr 2 + r 2 ^ 2 and so 

K = j-?yQ 2 Y~ 2 . We have 

{ 4e 2 l _ 1 

r(i+B a ) R — ' 

r& |fc-Z|=l 

else 

which has the spectrum 

with p = 1 . . . n. The nonzero eigenvalue smallest in magnitude is so for large n and small 0r _1 , approximately 

26 2 /2tt^ 2 

soa = ^r(^) 2 and 



r /ajr> 
V n 
o \ 2 

2n 



\fl + E 2 . 

Moving the numbers into constants k\ and £2, and dropping the 1 in 1 + /3 (fine for large n), we have 



At < kiVl + E^Y- 1 ! 



4 



AT rc i < e 2 r~ 2 n 2 . 

We again estimate \\f\\, now each column and row of v has 2 or 4 nonzero entries and so 

VlQ < IMI < V2® 

for some scaling and dimension independent constants v\ and v^. Then Theorem [5] gives the bound 

Ari.rd < ^v 2 n 2 QY- 1 . 

This time the condition 

M\<\\\b- l V 

does not break down for large dimensions, so the bound holds for all n when and Y are kept constant. 
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